JOURNAL OF RESEARCH of the National Bureau of Standards- A. Physics and Chemistry 
Vol. 73A, No. 1, January-February 1969 

Morphological Stability of a Cylinder 

Sam R. Coriell and Stephen C. Hardy 

Institute for Materials Research, National Bureau of Standards, Washington, D.C. 20234 

(August 27, 1968) 

The stability of the shape of a solid cylinder crystallizing in a supercooled liquid is treated. The 
effects of solute diffusion, slightly anisotropic surface tension and interface kinetics are included. The 
resulting stability equations are applied to the specific case of ice cylinders. 
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1. Introduction 

The stability of the shape of a solid growing by dif- 
fusion or heat flow was first studied by Mullins and 
Sekerka [l] 1 , who determined the stability criteria 
for a sphere. The cylindrical geometry has been 
treated by Coriell and Parker [2] and Kotler and 
Tiller [3]. These studies assumed that the interface 
properties were isotropic. Cahn [4] has taken account 
of slightly anisotropic surface free energy for the 
sphere. Recently, we have made experimental studies 
[5, 6] of the morphological stability of ice cylinders. 
In connection with this experimental program, it 
appears desirable to work out the stability of a cylin- 
drical shape taking account of slightly anisotropic 
surface tension and interface kinetics. In this paper 
we treat the case of a cylinder crystallizing from a 
binary melt; it is assumed that the surface tension 
and kinetic coefficient are slightly anisotropic. For 
the case of isotropic interface properties, the problem 
reduces to the case previously treated by Kotler and 
Tiller [3]. 

Although our calculation is general, we are par- 
ticularly interested in the case of an ice cylinder 
growing from pure water and from water with im- 
purities added. We present some specific calculations 
for these cases. 

2. Formulation and Calculation 

We wish to solve Laplace's equation [1, 2] for the 
temperature T and concentration C of impurity. The 
subscripts S and L denote solid and liquid, respec- 
tively, while a subscript / denotes that the quantity 



is to be evaluated at the interface. We assume a slightly 
perturbed cylinder of shape [2] 



r = /? + 8e i/c V 277 ^ x , 



(1) 



where r, <£, and z are the usual cylindrical coordi- 
nates, R is the radius of the unperturbed cylinder, 
8 is the amplitude of the perturbation (8/R < 1), and 
k and k determine the shape of the perturbation. 
Neglecting diffusion in the solid, the equations to be 
solved are 



with the auxiliary equations 

T L (R b ) = T b 
T s (0) -finite 

Tu = T S i 
Cl(Rc) =Cb 



ks/dTs 
U \dr n 



k(§Ik 



»- {2 /*****} {Te-Tu} 



= -D L (dC h 
' C u (l-j) \dr„ 



1 Figures in brackets indicate the literature references at the end of this paper. 



T e = T M - T M K £ !>'** - mCu 



(2) 

(3a) 
(3b) 
(3c) 
(3d) 

(3e) 
(30 

(3g) 

(3h) 
(3i) 
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and 



H, = kJ' k {k z )lh{kz), 



1 Fi 

K= R^~R* ( k2 + k2 z- l ) eik<Pei2nzlX - ( 3 J) 

In the above equations Rb is the bath radius, Tb (a 
constant) is the temperature at Rb, R c is the radius 
at which the concentration of impurity is C& (a con- 
stant), v is the velocity of the interface, k s and k L are 
the thermal conductivities of the solid and liquid, 
respectively, L v is the latent heat per unit volume of 

the solid, the f V jJL k e ik< ^ I is the slightly anisotropic 

linear kinetic coefficient, 2 T e is the equilibrium tem- 
perature of the interface, Dl is the solute diffusion 
coefficient, j is the partition coefficient (the ratio of the 
equilibrium concentration of solute on the solid side 
of the interface to that on the liquid side of the inter- 
face), R and 8 are the rate of growth of the unper- 
turbed cylinder and the perturbation, respectively, Tm 
is the melting point of a flat interface, K is the curva- 
ture of the perturbed cylinder, the f^r^e^l is the 

slightly anisotropic capillary constant, m is the freez- 
ing point lowering constant, and k z =27rRIX. The solu- 
tion of Laplace's equation in cylindrical coordinates 
is of the form C,-hC 2 In r+ e**e*«*l x \CJ k (k t rlR) 
+ C4Kk(k z r/R)] where the C's are constants and h 
and Kk are modified Bessel functions. A solution of 
this form is written for 7l, r s , and Cl and the con- 
stants are determined by the boundary conditions. 
Since we are treating small perturbations and slightly 
anisotropic surface tension and interface kinetics, 
we neglect all terms that are greater than first order 
in 8, fik, and V k (k # 0), i.e., we omit terms containing 
8 2 , 8/jik, STu, TkHk, etc., (k # 0). Although the calcu- 
lation is tedious, it is straightforward and we will not 
reproduce the details here. We define Ab = \n (RblR), 
A c = ln (R c /R), 

bT=TM-{T u YolR)-T b , A6=AT-mC b 
p A = (kJRLvfM,) +A b , f = [0.-j)hAcV[DL,] 

A,= [K k {k z R b IR)lh{k z R b IR)], 

D x =[K k (kzR c IR)lh(k z R c IR)l 

jA = -k z {[K k (k z )-AJ k (k z )]l[K k (k z )-AJ k (k z )]}, 

J» = -k z {[K' k {k z )-DJ' k {k z )]l[K k {k z )-DJ k {k z )]}, 



2 The kinetic coefficient and the capillary constant can be written in this form since we 
are only considering first order terms in the perturbation amplitude and in the anisotropy. 
In general, for example, the kinetic coefficient /x can be expanded as 



fl = ^ [ dk z ix k {k z )e ik *'e ik z { *' !R 



where <P' and z' give the orientation of the surface. For a slightly perturbed cylinder, the 
orientational coordinates z and <$>' are related to the space coordinates z and <t> by z'=0+8a t 
and <£>' = <$> + 8a 2 , where a t and a% are functions of <t> and z. Substituting these expressions 
in the expansion of /x, expanding the exponentials, neglecting terms of the order fi. k (k z )8, 
and performing the integral over k z we obtain fi = 1fji k e ik< ', where fji k = ^dk z fx k (k z ). A sim- 
ilar procedure gives Y = V Y k e ik4> ; here Y is (1/Lr)(y + rf 2 y/d3> 2 ) where y is the surface tension. 



where the prime indicates the derivative with respect 
to k z . The results of the calculation are 



R=(k I JRL v A b )(DT), 
where (DT) is given by 



(DT)=A b 



_1_ AT 

2^2p A 



J_ AT 

2^2i3 A 



(4a) 
(4b) 



For the ratio (S/8)/ (#//?) we obtain 



l={ 1+ £ + v}> 



-i+=^Cfo-D 



[jA-T-ikslk^H^Ab lkRiuLk 
(DT) [ 8tf 

T M T k T M T«(k 2 + k 2 z -l) 



8 R 

where a = (k L /L v R ) (J a + k s Hilk L ) 
and 



(4c) 



kiAc 
LyRCbt; 



1 



{(DT) 



Ab 



Jd- 



HDT) 



A c A b 



In the above equations F k and /^(A:#0) are zero if 
k z ^ 0, i.e., the anisotropic terms only effect O-type 
perturbations and have no effect on perturbations 
along the axis of the cylinder. This follows, for example, 
since there are no terms of the form T k , kz e ik ^e i27TzfX in 
eq (3i) (see footnote 2). 

Thus one of the important results of this calculation 
is that the growth rate of z-perturbations is independent 
of slightly anisotropic surface tension and kinetic 
coefficient. 

Discussion 

The general results of our stability analysis are given 
by eqs (4a-c). In this section, we apply these equations 
to the special case of ice cylinders (oriented with the 
c-axis of ice parallel to the cylinder axis) growing in 
slightly supercooled water {AT = 0.1 °C). We also 
consider the case in which impurities have been added 
to the water. We are interested in the case where R 
lies between 0.04 cm and 0.25 cm; the bath radius 
Rb in most of the experimental work was 0.875 cm. 
Under the above conditions we can simplify the equa- 
tions by neglecting interface kinetics, i.e., we take 

l+-sl. This implies that (kJRLv) (J U + H,k s lk L ) 
Mo 

< Mo; typically (J A + #jfe/fa,) < 100 and the above 

approximation is valid if /x () > .05 cm deg -1 sec -1 . 
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Although the interface kinetic coefficient for ice grow- 
ing perpendicular to the c-axis is not known, a lower 
bound can be obtained from bath undercooling meas- 
urements. For example, Lindenmeyer et al. [7], ob- 
served a growth velocity of 1 cm sec -1 at a bath under- 
cooling of 2 °C. Assuming a linear law, this gives 
/jlo > .5 cm sec -1 deg -1 [8]. Since fio may be orders 
of magnitude greater than this, it appears valid to 
neglect interface kinetics. It follows that (3a = A^. 

It is almost always (and certainly under the experi- 
mental conditions of interest) valid to neglect the A\ 
and D\ terms in the definition of J a and Jd. Defining 
H K = — k z K' k (k z )l(K k (k z ), we then have J A =H k =Jd. 
With this approximation and the neglect of interface 
kinetics, eqs (4b) and (4c) can be written 



(DT) = i( [ (4M) + AT] - | { [At/f) + AT] 2 

-4Afc4»/£}»/*|) 



(5a) 



and 



8R -H 
SR-"*' 



t [H K + (kalkAHAAt 

{\ + majv)DT 



X L 8,4 ' 8 R J' (bb) 

For pure water, there is further simplification, viz, 
taking m=0 and f— *0, 



(DT) = AT 



(6a) 



and 



8R-h i i [MjLtMklMlMi 
8R' Mk 1+ DT 



RRn k T M r k T M U(k 2 + kl 
8fi 2 8 R 



ill 



(6b) 



Recalling that for a k z ^ perturbation /i* = and 
I\ = (k 7^ 0), eq (6b) becomes 



8R_ 
8R 



H K -l 



T M r (k 2 +ki 



■l)A b \H K - 
R(DT) 



(k s /k L )HA 



(6c) 



This equation has been used in references 5 and 6 to 
analyze the experimental data. 

We now give a brief discussion of the choice of the 
bath boundary conditions, i.e., the choice of Ri, and 
R c ; we also discuss the use of Laplace's equation. For 
an infinite bath, the choice of these parameters has 
been previously discussed [2]. For an infinite bath 
Rt, = Rx where Rx = /?/(1.33X) and A satisfies 

k 2 e k2 Ei(-k 2 )+S = 0. 



In the above equation Ei is the exponential integral 
function and for heat flow S = St = C v (T/ — Ti ) )IL v and 
for impurity diffusion [9] S = S C = (Cu"~C Q )ICuO —j). 
In these equations T^ and Cf, are the temperature and 
concentration at infinity, respectively, and C v is the 
specific heat per unit volume of the liquid. In order for 
Laplace's equation to be valid, it is necessary that 

St < 1 and S c ^ 1. Since St= (^t) A7 1 , it is clear that 

for A7< 1, St<^ 1. Using the solution of the diffusion 
equation, we may rewrite S c as S c = {A c IA b )(kJDL v )DT. 
Taking fe = 1.33(10- 3 ) cal cm" 1 deg" 1 sec" 1 , L,= 73.4 
cal cm -3 , D= 10~ 5 cm 2 sec -1 , and (Ac/At) = 1, we have 
S c = 1.8 (DT). For DT =0.1, S c = 0.2, and Laplace's 
equation is a reasonable approximation. For larger 
{DT) and for smaller diffusion constants, there may 
well be deviations from Laplace's equation. 

In the actual experiments the temperature is 
maintained at Tt> at some radius, say R a . If R a < Rx, 
it seems reasonable to take Ri, = R a and this has been 
done [5, 6]. If R a > Rx, it is better to take R h = Rx. If 
this were not done, i.e., if instead we let Ri, = R„ when 
R a > Rx, then from eq (4a) we would predict that 
R is larger when TiX QO ) = T b than when TiXR a ) — Ti). 
This is obviously wrong and hence for R„ > Rx we take 
Rb = Rx. When R a < Rx, there is a certain error in 
taking R b = R a , but this error lies within the experi- 
mental error in measuring AT and R. The choice of 
R c is even more complicated since the average con- 
centration in the liquid changes as solute is rejected 
from the ice. Fortunately for many cases the results 
are not very sensitive to the choice of R c and for 
calculational purposes we take R c = Ra- The proper 
choice of R c and Rb can be studied experimentally 
by comparing the observed growth rate ft of the 
unperturbed cylinder with the theoretical R. 

We present some calculations of 8 as a function of 
R for the special case where (8RI8R) is given by eq 
(6c). Following Cahn [4] we use the relationship 
(8RI8k) = (RI8)(d8ldR). Denoting the right hand 
side of eq (6c) by f(R), we have 



d In 8f(R) 



dR 



R 



Integrating yields 



In (81 '& 



Jfto X 



x) 



dx, 



(7a) 



(7b) 



where 8 is the value of 8 at R = R {) . The integral has 
been evaluated numerically for various values of the 
parameters occurring in f(R). In many cases over 
small ranges of R, In 8 is to a good approximation 
linear with R. This is illustrated in table 1 which 
gives (d\n8/dR) as a function of R. From table 1, 
it is seen that (din 8/dR) changes very slowly with R. 
Since in the experimental measurements of 8 versus/?, 
R varies by about 0.02 to 0.03 cm and R is greater 
than 0.120 cm, it is to be expected that plots of the 
experimental data in the form In 8 versus R will 
appear linear. 
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Of interest is the wavelength dependence of the 
8 versus R curves. This is shown in table 2 for three 
different wavelengths. We can consider 8/80 as an 
amplification factor, that is, it gives the relative 
magnitude of 8 at R as compared to 8 at R . As seen 
from the table, these amplification factors are very 
large. Experimentally values of 8 of 3(10~ 3 ) cm have 
been observed at /?= 0.175 cm. Extrapolating to 
/? = 0.05 cm gives a value of 8 = 10~ 7 cm as the size 
of the initial perturbation. It is clear from table 2 
that the amplification factor is a function of wavelength. 
For small R, perturbations with X= 0.045 cm are ampli- 
fied slightly more than perturbations corresponding 
to the other wavelengths of the table. For larger /?, 
however, the amplification is greatest for A. = 0.055. 
Calculations similar to these will be useful in attempt- 
ing to predict the wavelength of the perturbation which 
appears on a growing cylinder. However, such a 
prediction requires some assumption about the 
initial distribution of the amplitudes of the perturba- 
tions of various wavelengths. 

We wish to discuss the effect of solute on the 
stability equations. It is interesting that for no aniso- 
tropy and zero surface tension, the addition of solute 
has no effect on (8/8)/ (/J//?), i.e., eqs (5b) and (6b) 
are identical. The addition of solute changes the 
unperturbed growth rate by changing DT. For a bath 
undercooling Ar=0.1, the value of (DT) for various 
values of mcb are given in table 3. 

Also of interest is the factor l + (ma/v) appearing 
in eq (5b). It can be shown that 

Mi/) = (kJDL P )(mc b )(l + k s /k L ) = 9 (mc b ). 

Thus for met, < 0.1 deg., \ + (malv) varies from 1 to 2. 
Since DT varies more rapidly than this, the main effect 
of adding solute is to decrease DT and thus to increase 
the magnitude of the last term in eq (5b). Thus for 
fixed A7 7 adding solute makes (8/8)/(/?//?) smaller, thus 
stabilizing the cylinder. On the other hand for fixed 
(DT), which corresponds to not changing the unper- 
turbed velocity /?, the addition of solute increases 
(&ld)l(R/R) and consequently makes the cylinder more 
unstable. The above statements assume that the 
To term is the dominant term inside the bracket of 
eq (5b) and that the addition of solute does not change 
the physical properties such as the surface tension. 
If solute is adsorbed at the interface, one expects a 
lowering of the surface tension. One method of study- 
ing such effects is to add a very small amount of solute, 
e.g., 10" 4 M NaCl in which case mc b = 4(10" 4 ). From the 
preceding calculations, it is clear that DT = A7 and 
1 + (majv) = 1 so that eqs (6a) and (6b) can be used. 
Thus any difference between experimental results 
for distilled water and water containing 10~ 4 M of 
impurities should probably be attributed to a change 
in the surface tension due to adsorption at the interface. 
In summary, we have analyzed the stability of a 
solid cylinder growing by heat flow into a binary melt 
and have taken account of any small anisotropy in 
the interface properties. In particular, we have pro- 
vided a theoretical framework for the experimental 
study of the stability of ice cylinders. 



Table 1. (d In 8/dR) as a function of R for \ = 0.055 cm; y = 0.018 
J/m 2 and AT '=0.1 deg. 



R(cm) 


(d\n8/dR) (cm- 1 ) 


0.100 


85.8 


.110 


84.3 


.120 


82.9 


.130 


81.7 


.140 


80.7 


.150 


79.7 


.160 


78.9 


.170 


78.1 


.180 


77.4 


.190 


76.8 


.200 


76.3 



Table 2. Effect of wavelength on 8/80 for y— 0.018 J/m 2 and 
AT =0.1 deg. 







(8/80) 




R(cm) 










X = 0.045 cm 


X = 0.055 cm 


X = 0.065 cm 


0.050 


1.00 


1.00 


1.00 


.075 


1.08(100 


1.07(100 


9.88 


.100 


9.85(100 


9.70(100 


8.16(100 


.125 


7.80(H) 2 ) 


7.91 (10 2 ) 


6.07(10 2 ) 


.150 


5.53(H) 3 ) 


5.99(10 3 ) 


4.22(10 3 ) 


.175 


3.59(H) 4 ) 


4.28(10 4 ) 


2.80(H) 4 ) 


.200 


2.17(H) 5 ) 


2.93(10 5 ) 


1.80(100 



Table 3. Value of effective undercooling (DT) as a function of added 
solute for bath undercooling AT = 0.1 °C (£/A b = 7.8 deg' 1 used in 
the calculation) 



mc b 


DT 


0.001 


0.0988 


.01 


.0881 


.02 


.0768 


.04 


.0556 


.06 


.0359 


.08 


.0174 
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